*! ranktest 1.2.03  24Apr2010
*! author mes, based on code by fk

program define ranktest, rclass sortpreserve
	version 9.2
	local lversion 01.2.03

	if substr("`1'",1,1)== "," {
		if "`2'"=="version" {
			di in ye "`lversion'"
			return local version `lversion'
			exit
		}
		else {
di as err "invalid syntax"
			exit 198
		}
	}

* If varlist 1 or varlist 2 have a single element, parentheses optional

	if substr("`1'",1,1)=="(" {
		GetVarlist `0'
		local y `s(varlist)'
		local K : word count `y'
		local 0 `"`s(rest)'"'
		sret clear
	}
	else {
		local y `1'
		local K 1
		mac shift 1
		local 0 `"`*'"'
	}

	if substr("`1'",1,1)=="(" {
		GetVarlist `0'
		local z `s(varlist)'
		local L : word count `z'
		local 0 `"`s(rest)'"'
		sret clear
	}
	else {
		local z `1'
		local K 1
		mac shift 1
* Need to reinsert comma before options (if any) for -syntax- command to work
		local 0 `", `*'"'
	}

* Option version ignored here if varlists were provided
	syntax [if] [in] [aw fw pw iw/] [, partial(varlist ts) fwl(varlist ts) /*
		*/ NOConstant wald ALLrank NULLrank FULLrank ROBust cluster(varlist) /*
		*/ BW(string) kernel(string) Tvar(varname) Ivar(varname) sw spsd version]

	local partial "`partial' `fwl'"

	marksample touse

	if "`noconstant'"=="" {
		tempvar one
		gen byte `one' = 1
		local partial "`partial' `one'"
	}

	if "`wald'"~="" {
		local LMWald "Wald"
	}
	else {
		local LMWald "LM"
	}
	
	local optct : word count `allrank' `nullrank' `fullrank'
	if `optct' > 1 {
di as err "Incompatible options: `allrank' `nullrank' `fullrank'"
		error 198
	}
	else if `optct' == 0 {
* Default
		local allrank "allrank"
	}

* Note that by tsrevar-ing here, subsequent disruption to the sort doesn't matter
* for TS operators.
	tsrevar `y'
	local vl1 `r(varlist)'
	tsrevar `z'
	local vl2 `r(varlist)'
	tsrevar `partial'
	local partial `r(varlist)'

	foreach vn of varlist `vl1' {
		tempvar tv
		qui gen double `tv' = .
		local tempvl1 "`tempvl1' `tv'"
	}
	foreach vn of varlist `vl2' {
		tempvar tv
		qui gen double `tv' = .
		local tempvl2 "`tempvl2' `tv'"
	}

* Stock-Watson and cluster imply robust.
	if "`sw'`cluster'" ~= "" {
		local robust "robust"
	}

	tempvar wvar
	if "`weight'" == "fweight" | "`weight'"=="aweight" {
		local wtexp `"[`weight'=`exp']"'
		gen double `wvar'=`exp'
	}
	if "`weight'" == "fweight" & "`kernel'" !="" {
		di in red "fweights not allowed (data are -tsset-)"
		exit 101
	}
	if "`weight'" == "fweight" & "`sw'" != "" {
		di in red "fweights currently not supported with -sw- option"
		exit 101
	}
	if "`weight'" == "iweight" {
		if "`robust'`cluster'`bw'" !="" {
			di in red "iweights not allowed with robust, cluster, AC or HAC"
			exit 101
		}
		else {
			local wtexp `"[`weight'=`exp']"'
			gen double `wvar'=`exp'
		}
	}
	if "`weight'" == "pweight" {
		local wtexp `"[aweight=`exp']"'
		gen double `wvar'=`exp'
		local robust "robust"
	}
	if "`weight'" == "" {
* If no weights, define neutral weight variable
		qui gen byte `wvar'=1
	}
	else if "`weight'"=="aweight" | "`weight'" == "pweight" {
		sum `wvar' if `touse', meanonly
		qui replace `wvar'=`wvar'*r(N)/r(sum)
	}

	markout `touse' `vl1' `vl2' `partial' `cluster', strok

* HAC estimation.
* If bw is omitted, default `bw' is empty string.
* If bw or kernel supplied, check/set `kernel'.
* Macro `kernel' is also used for indicating HAC in use.
	if "`bw'" != "" | "`kernel'" != "" {
* Need tvar for markout with time-series stuff
* Data must be tsset for time-series operators in code to work
* User-supplied tvar checked if consistent with tsset
		capture tsset
		if "`r(timevar)'" == "" {
di as err "must tsset data and specify timevar"
			exit 5
		}
		if "`tvar'" == "" {
			local tvar "`r(timevar)'"
		}
		else if "`tvar'"!="`r(timevar)'" {
di as err "invalid tvar() option - data already -tsset-"
			exit 5
		}
* If no panel data, ivar will still be empty
		if "`ivar'" == "" {
			local ivar "`r(panelvar)'"
		}
		else if "`ivar'"!="`r(panelvar)'" {
di as err "invalid ivar() option - data already -tsset-"
			exit 5
		}
		local tdelta `r(tdelta)'
		tsreport, panel
		if `r(N_gaps)' != 0 {
di in gr "Warning: time variable " in ye "`tvar'" in gr " has " /*
	*/ in ye "`r(N_gaps)'" in gr " gap(s) in relevant range"
		}
		if "`bw'" == "" {
di as err "bandwidth option bw() required for HAC-robust estimation"
			exit 102
		}
		local bw real("`bw'")
* Check it's a valid bandwidth; allow non-integer bandwidth
//		if   `bw' != int(`bw') | /*
		if   `bw' == .  | `bw' <= 0 {
di as err "invalid bandwidth in option bw() - must be integer > 0"
			exit 198
		}
* Convert bw macro to simple integer
		local bw=`bw'

* Check it's a valid kernel
		local validkernel 0
		if lower(substr("`kernel'", 1, 3)) == "bar" | "`kernel'" == "" {
* Default kernel
			local kernel "Bartlett"
			local window "lag"
			local validkernel 1
			if `bw'==1 {
di in ye "Note: kernel=Bartlett and bw=1 implies zero lags used."
di in ye "      Test statistics are not autocorrelation-consistent."
			}
		}
		if lower(substr("`kernel'", 1, 3)) == "par" {
			local kernel "Parzen"
			local window "lag"
			local validkernel 1
			if `bw'==1 {
di in ye "Note: kernel=Parzen and bw=1 implies zero lags used."
di in ye "      Test statistics are not autocorrelation-consistent."
			}
		}
		if lower(substr("`kernel'", 1, 3)) == "tru" {
			local kernel "Truncated"
			local window "lag"
			local validkernel 1
		}
		if lower(substr("`kernel'", 1, 9)) == "tukey-han" | lower("`kernel'") == "thann" {
			local kernel "Tukey-Hanning"
			local window "lag"
			local validkernel 1
			if `bw'==1 {
di in ye "Note: kernel=Tukey-Hanning and bw=1 implies zero lags."
di in ye "      Test statistics are not autocorrelation-consistent."
			}
		}
		if lower(substr("`kernel'", 1, 9)) == "tukey-ham" | lower("`kernel'") == "thamm" {
			local kernel "Tukey-Hamming"
			local window "lag"
			local validkernel 1
			if `bw'==1 {
di in ye "Note: kernel=Tukey-Hamming and bw=1 implies zero lags."
di in ye "      Test statistics are not autocorrelation-consistent."
			}
		}
		if lower(substr("`kernel'", 1, 3)) == "qua" | lower("`kernel'") == "qs" {
			local kernel "Quadratic spectral"
			local window "spectral"
			local validkernel 1
		}
		if lower(substr("`kernel'", 1, 3)) == "dan" {
			local kernel "Daniell"
			local window "spectral"
			local validkernel 1
		}
		if lower(substr("`kernel'", 1, 3)) == "ten" {
			local kernel "Tent"
			local window "spectral"
			local validkernel 1
		}
		if ~`validkernel' {
			di in red "invalid kernel"
			exit 198
		}
	}
	else {
		local bw 1
	}

* tdelta missing if version 9 or if not tsset			
	if "`tdelta'"=="" {
		local tdelta=1
	}

	if "`sw'"~="" {
		capture xtset
		if "`ivar'" == "" {
			local ivar "`r(panelvar)'"
		}
		else if "`ivar'"!="`r(panelvar)'" {
di as err "invalid ivar() option - data already tsset or xtset"
			exit 5
		}
* Exit with error if ivar is neither supplied nor tsset nor xtset
		if "`ivar'"=="" {
di as err "Must -xtset- or -tsset- data or specify -ivar- with -sw- option"
			exit 198
		}
		qui describe, short varlist
		local sortlist "`r(sortlist)'"
		tokenize `sortlist'
		if "`ivar'"~="`1'" {
di as err "Error - dataset must be sorted on panel var with -sw- option"
			exit 198
		}
	}

* Create variable used for getting lags etc. in Mata
	tempvar tindex
	qui gen `tindex'=1 if `touse'
	qui replace `tindex'=sum(`tindex') if `touse'

********** CLUSTER SETUP **********************************************

* Mata code requires data are sorted on (1) the first var cluster if there
* is only one cluster var; (2) on the 3rd and then 1st if two-way clustering,
* unless (3) two-way clustering is combined with kernel option, in which case
* the data are tsset and sorted on panel id (first cluster variable) and time
* id (second cluster variable).
* Second cluster var is optional and requires an identifier numbered 1..N_clust2,
* unless combined with kernel option, in which case it's the time variable.
* Third cluster var is the intersection of 1 and 2, unless combined with kernel
* opt, in which case it's unnecessary.
* Sorting on "cluster3 cluster1" means that in Mata, panelsetup works for
* both, since cluster1 nests cluster3.
* Note that it is possible to cluster on time but not panel, in which case
* cluster1 is time, cluster2 is empty and data are sorted on panel-time.
* Note also that if no kernel-robust, sorting will disrupt any tsset-ing,
* but data are tsrevar-ed earlier to avoid any problems.
	if "`cluster'"!="" {
		local clopt "cluster(`cluster')"
		tokenize `cluster'
		local cluster1 "`1'"
		local cluster2 "`2'"
		if "`kernel'"~="" {
* kernel requires either that cluster1 is time var and cluster2 is empty
* or that cluster1 is panel var and cluster2 is time var.
* Either way, data must be tsset and sorted for panel data.
			if "`cluster2'"~="" {
* Allow backwards order
				if "`cluster1'"=="`tvar'" & "`cluster2'"=="`ivar'" {
					local cluster1 "`2'"
					local cluster2 "`1'"
				}
				if "`cluster1'"~="`ivar'" | "`cluster2'"~="`tvar'" {
di as err "Error: cluster kernel-robust requires clustering on tsset panel & time vars."
di as err "       tsset panel var=`ivar'; tsset time var=`tvar'; cluster vars=`cluster1',`cluster2'"
					exit 198
				}
			}
			else {
				if "`cluster1'"~="`tvar'" {
di as err "Error: cluster kernel-robust requires clustering on tsset time variable."
di as err "       tsset time var=`tvar'; cluster var=`cluster1'"
					exit 198
				}
			}
		}
* Simple way to get quick count of 1st cluster variable without disrupting sort
* clusterid1 is numbered 1.._Nclust1.
		tempvar clusterid1
		qui egen `clusterid1'=group(`cluster1') if `touse'
		sum `clusterid1' if `touse', meanonly
		if "`cluster2'"=="" {
			local N_clust=r(max)
			local N_clust1=.
			local N_clust2=.
			if "`kernel'"=="" {
* Single level of clustering and no kernel-robust, so sort on single cluster var.
* kernel-robust already sorted via tsset.
				sort `cluster1'
			}
		}
		else {
			local N_clust1=r(max)
			if "`kernel'"=="" {
				tempvar clusterid2 clusterid3
* New cluster id vars are numbered 1..N_clust2 and 1..N_clust3
				qui egen `clusterid2'=group(`cluster2') if `touse'
				qui egen `clusterid3'=group(`cluster1' `cluster2') if `touse'
* Two levels of clustering and no kernel-robust, so sort on cluster3/nested in/cluster1
* kernel-robust already sorted via tsset.
				sort `clusterid3' `cluster1'
			}
			else {
				local clusterid2 `cluster2'
			}
			sum `clusterid2' if `touse', meanonly
			local N_clust2=r(max)
			local N_clust=min(`N_clust1',`N_clust2')
		}
	}

************************************************************************************************


* Note that bw is passed as a value, not as a string
	mata: rkstat(	"`vl1'",			/*
				*/	"`vl2'",			/*
				*/	"`partial'",		/*
				*/	"`wvar'",			/*
				*/	"`weight'",			/*
				*/	"`touse'",			/*
				*/	"`LMWald'",			/*
				*/	"`allrank'",		/*
				*/	"`nullrank'",		/*
				*/	"`fullrank'",		/*
				*/	"`robust'",			/*
				*/	"`clusterid1'",		/*
				*/	"`clusterid2'",		/*
				*/	"`clusterid3'",		/*
				*/	`bw',				/*
				*/	"`tvar'",			/*
				*/	"`ivar'",			/*
				*/	"`tindex'",			/*
				*/	`tdelta',			/*
				*/	"`kernel'",			/*
				*/	"`sw'",				/*
				*/	"`spsd'",			/*
				*/	"`window'",			/*
				*/	"`tempvl1'",		/*
				*/	"`tempvl2'")

	tempname rkmatrix chi2 df df_r p rank ccorr eval
	mat `rkmatrix'=r(rkmatrix)
	mat `ccorr'=r(ccorr)
	mat `eval'=r(eval)
	mat colnames `rkmatrix' = "rk" "df" "p" "rank" "eval" "ccorr"
	
di
di "Kleibergen-Paap rk `LMWald' test of rank of matrix"
	if "`robust'"~="" & "`kernel'"~= "" & "`cluster'"=="" {
di "  Test statistic robust to heteroskedasticity and autocorrelation"
di "  Kernel: `kernel'   Bandwidth: `bw'"
	}
	else if "`kernel'"~="" & "`cluster'"=="" {
di "  Test statistic robust to autocorrelation"
di "  Kernel: `kernel'   Bandwidth: `bw'"
	}
	else if "`cluster'"~="" {
di "  Test statistic robust to heteroskedasticity and clustering on `cluster'"
		if "`kernel'"~="" {
di "  and kernel-robust to common correlated disturbances"
di "  Kernel: `kernel'   Bandwidth: `bw'"
		}
	}
	else if "`robust'"~="" {
di "  Test statistic robust to heteroskedasticity"
	}
	else if "`LMWald'"=="LM" {
di "  Test assumes homoskedasticity (Anderson canonical correlations test)"
	}
	else {
di "  Test assumes homoskedasticity (Cragg-Donald test)"
	}
		
	local numtests = rowsof(`rkmatrix')
	forvalues i=1(1)`numtests' {
di "Test of rank=" %3.0f `rkmatrix'[`i',4] "  rk=" %8.2f `rkmatrix'[`i',1] /*
	*/	"  Chi-sq(" %3.0f `rkmatrix'[`i',2] ") pvalue=" %8.6f `rkmatrix'[`i',3]
	}
	scalar `chi2' = `rkmatrix'[`numtests',1]
	scalar `p' = `rkmatrix'[`numtests',3]
	scalar `df' = `rkmatrix'[`numtests',2]
	scalar `rank' = `rkmatrix'[`numtests',4]
	local N `r(N)'
	return scalar df = `df'
	return scalar chi2 = `chi2'
	return scalar p = `p'
	return scalar rank = `rank'
	if "`cluster'"~="" {
		return scalar N_clust = `N_clust'
	}
	if "`cluster2'"~="" {
		return scalar N_clust1 = `N_clust1'
		return scalar N_clust2 = `N_clust2'
	}
	return scalar N = `N'
	return matrix rkmatrix `rkmatrix'
	return matrix ccorr `ccorr'
	return matrix eval `eval'
	
	tempname S V Omega
	if `K' > 1 {
		foreach en of local y {
* Remove "." from equation name
			local en1 : subinstr local en "." "_", all
			foreach vn of local z {
				local cn "`cn' `en1':`vn'"
			}
		}
	}
	else {
		foreach vn of local z {
		local cn "`cn' `vn'"
		}
	}
	mat `V'=r(V)
	matrix colnames `V' = `cn'
	matrix rownames `V' = `cn'
	return matrix V `V'
	mat `S'=r(S)
	matrix colnames `S' = `cn'
	matrix rownames `S' = `cn'
	return matrix S `S'	
end

* Adopted from -canon-
program define GetVarlist, sclass 
	sret clear
	gettoken open 0 : 0, parse("(") 
	if `"`open'"' != "(" {
		error 198
	}
	gettoken next 0 : 0, parse(")")
	while `"`next'"' != ")" {
		if `"`next'"'=="" { 
			error 198
		}
		local list `list'`next'
		gettoken next 0 : 0, parse(")")
	}
	sret local rest `"`0'"'
	tokenize `list'
	local 0 `*'
	sret local varlist "`0'"
end

version 9.2
mata:
void rkstat(	string scalar vl1,
				string scalar vl2,
				string scalar partial,
				string scalar wvarname,
				string scalar weight,
				string scalar touse,
				string scalar LMWald,
				string scalar allrank,
				string scalar nullrank,
				string scalar fullrank,
				string scalar robust,
				string scalar clustvarname,
				string scalar clustvarname2,
				string scalar clustvarname3,
				bw,
				string scalar tvarname,
				string scalar ivarname,
				string scalar tindexname,
				tdelta,
				string scalar kernel,
				string scalar sw,
				string scalar spsd,
				string scalar window,
				string scalar tempvl1,
				string scalar tempvl2)
{

// tempx, tempy and tempz are the Stata names of temporary variables that will be changed by rkstat
	if (partial~="") {
		tempx=tokens(partial)
	}
	tempy=tokens(tempvl1)
	tempz=tokens(tempvl2)

	st_view(y=.,.,tokens(vl1),touse)
	st_view(z=.,.,tokens(vl2),touse)
	st_view(yhat=.,.,tempy,touse)
	st_view(zhat=.,.,tempz,touse)
	st_view(mtouse=.,.,tokens(touse),touse)
	st_view(wvar=.,.,tokens(wvarname),touse)
	noweight=(st_vartype(wvarname)=="byte")

// Effective number of observations is sum of weight variable.
// If no weights, then wvar is ones.
	N=colsum(wvar)

	if (clustvarname~="") {
		st_view(clustvar,.,clustvarname,touse)
		info = panelsetup(clustvar, 1)
		N_clust=rows(info)
		if (clustvarname2~="") {
			st_view(clustvar2, ., clustvarname2, touse)
			if (kernel=="") {
				st_view(clustvar3, ., clustvarname3, touse) // needed only if not panel tsset
			}
		}
	}

	if (kernel~="") {
		st_view(t, ., tvarname, touse)
		T=max(t)-min(t)+1
	}

// Partial out the X variables
	if (partial~="") {
		st_view(x=.,.,tempx,touse)
		xx = quadcross(x, wvar, x)
		xy = quadcross(x, wvar, y)
		xz = quadcross(x, wvar, z)
// y=XB => X'y=X'XB
		by = qrsolve(xx,xy)
		bz = qrsolve(xx,xz)
		yhat[.,.] = y-x*by
		zhat[.,.] = z-x*bz
	}
	else {
		yhat[.,.] = y
		zhat[.,.] = z
	}

	K=cols(y)
	L=cols(z)

	zhzh = quadcross(zhat, wvar, zhat)
	zhyh = quadcross(zhat, wvar, yhat)
	yhyh = quadcross(yhat, wvar, yhat)

	pihat = qrsolve(zhzh,zhyh)
// rzhat is F in paper (p. 103)
// iryhat is G in paper (p. 103)
	ryhat=cholesky(yhyh)
	rzhat=cholesky(zhzh)
	iryhat=luinv(ryhat')
	irzhat=luinv(rzhat')
	that=rzhat'*pihat*iryhat

// cc is canonical correlations.  Squared cc is eigenvalues.
	fullsvd(that, ut, cc, vt)
	vt=vt'
	vecth=vec(that)
	ev = cc:^2
// S matrix in paper (p. 100).  Not used in code below.
//	smat=fullsdiag(cc, rows(that)-cols(that))

	if (abs(1-cc[1,1])<1e-10) {
printf("\n{text:Warning: collinearities detected between (varlist1) and (varlist2)}\n")
	}
	if ((missing(ryhat)>0) | (missing(iryhat)>0) | (missing(rzhat)>0) | (missing(irzhat)>0)) {
printf("\n{error:Error: non-positive-definite matrix. May be caused by collinearities.}\n")
		exit(error(3351))
	}

// If Wald, yhat is residuals
	if (LMWald=="Wald") {
		yhat[.,.]=yhat-zhat*pihat
		yhyh = quadcross(yhat, wvar, yhat)
	}

// Covariance matrices
// vhat is W in paper (eqn below equation 17, p. 103)
// shat is V in paper (eqn below eqn 15, p. 103)
// shat * 1/N is same as estimated S matrix of orthog conditions as e.g. saved by ivreg2
//    (including cluster case)

// Start calculation of shat
	if ((robust=="") & (clustvarname=="")) {
// Block for homoskedastic and AC.
		sigma2=yhyh/N
		shat=sigma2#zhzh

		if (kernel~="") {
			if (window=="spectral") {
				TAU=T/tdelta-1
			}
			else {
				TAU=bw
			}
			tnow=st_data(., tindexname)
			for (tau=1; tau<=TAU; tau++) {
				kw = m_calckw(tau, bw, kernel)
				if (kw~=0) {						// zero weight possible with some kernels
													// save an unnecessary loop if kw=0
													// remember, kw<0 possible with some kernels!
					lstau = "L"+strofreal(tau)
					tlag=st_data(., lstau+"."+tindexname)
					tmatrix = tnow, tlag
					svar=(tnow:<.):*(tlag:<.)		// multiply column vectors of 1s and 0s
					tmatrix=select(tmatrix,svar)	// to get intersection, and replace tmatrix
// col 1 of tmatrix has row numbers of all rows of data with this time period that have a corresponding lag
// col 2 of tmatrix has row numbers of all rows of data with lag tau that have a corresponding ob this time period.
// if no lags exist, tmatrix has zero rows.  Code below still works (ghat+ghat'=0) but best to catch anyway
					if (rows(tmatrix)>0) {
// Should never happen that fweights or iweights make it here,
// but if they did the next line would be sqrt(wvari)*sqrt(wvari1)*wf
						wv=wvar[tmatrix[.,1]]:*wvar[tmatrix[.,2]]	// inner weighting matrix for quadcross
						sigmahat = quadcross(yhat[tmatrix[.,1],.], wv, yhat[tmatrix[.,2],.]) / N
						zzhat    = quadcross(zhat[tmatrix[.,1],.], wv, zhat[tmatrix[.,2],.])
						ghat = sigmahat#zzhat
						shat=shat+kw*(ghat+ghat')
					}
				}	// end non-zero kernel weight block
			}	// end tau loop
		}	// end kernel code
	}  // end homoskedastic and AC block

// Block for robust HC and HAC but not Stock-Watson and single clustering.
// Need to enter for double-clustering if one cluster is time.
	if ( (robust~="") & (sw=="") & ((clustvarname=="") | ((clustvarname2~="") & (kernel~="")))  ) {
// Block for HC and HAC
		shat=J(L*K,L*K,0)
		for (i=1; i<=rows(yhat); i++) {
			yzi=yhat[i,.]#zhat[i,.]
			if ((weight=="fweight") | (weight=="iweight")) {
// wvar is a column vector
				shat=shat+quadcross(yzi,yzi)*wvar[i]
			}
			else {
				shat=shat+quadcross(yzi,yzi)*(wvar[i])^2
			}
		}
		if (kernel~="") {
// Spectral windows require looping through all T-1 autocovariances
			if (window=="spectral") {
				TAU=T/tdelta-1
			}
			else {
				TAU=bw
			}
			tnow=st_data(., tindexname)
			for (tau=1; tau<=TAU; tau++) {
				ghat=J(L*K,L*K,0)
				kw = m_calckw(tau, bw, kernel)
				if (kw~=0) {						// zero weight possible with some kernels
													// save an unnecessary loop if kw=0
													// remember, kw<0 possible with some kernels!
					lstau = "L"+strofreal(tau)
					tlag=st_data(., lstau+"."+tindexname)
					tmatrix = tnow, tlag
					svar=(tnow:<.):*(tlag:<.)		// multiply column vectors of 1s and 0s
					tmatrix=select(tmatrix,svar)	// to get intersection, and replace tmatrix
// if tmatrix has zero rows, loop not entered
					for (i=1; i<=rows(tmatrix); i++) {
						wvari =wvar[tmatrix[i,1]]
						wvari1=wvar[tmatrix[i,2]]
						yi    =yhat[tmatrix[i,1],.]
						yi1   =yhat[tmatrix[i,2],.]
						zi    =zhat[tmatrix[i,1],.]
						zi1   =zhat[tmatrix[i,2],.]
						yzi =yi#zi
						yzi1=yi1#zi1
// Should never happen that fweights or iweights make it here, but if they did
// the next line would be ghat=ghat+yzi'*yzi1*sqrt(wvari)*sqrt(wvari1)
						ghat=ghat+quadcross(yzi,yzi1)*wvari*wvari1
					}
					shat=shat+kw*(ghat+ghat')
				}
			}
		}
	} // end HC and HAC block
	
	if (clustvarname~="") {
// Block for cluster-robust including double-clustering
// 2-level clustering: S = S(level 1) + S(level 2) - S(level 3 = intersection of levels 1 & 2)

// Prepare shat3 if 2-level clustering
		if (clustvarname2~="") {
			if (kernel~="") {		// second cluster variable is time
									// shat3 was already calculated above as shat
				shat3=shat
			}
			else {					// calculate shat3
									// data were sorted on clustvar3-clustvar1 so
									// clustvar3 is nested in clustvar1 and Mata panel functions
									// work for both.
				info3 = panelsetup(clustvar3, 1)
				if (rows(info3)==rows(yhat)) {		// intersection of levels 1 & 2 are all single obs
													// so no need to loop through row by row
					shat3=J(L*K,L*K,0)
					for (i=1; i<=rows(yhat); i++) {
						yzi=yhat[i,.]#zhat[i,.]
						shat3=shat3+quadcross(yzi,yzi)*(wvar[i])^2
					}
				}
				else {								// intersection of levels 1 & 2 includes some groups of obs
					N_clust3=rows(info3)
					shat3=J(L*K,L*K,0)
					for (i=1; i<=N_clust3; i++) {
						yz=J(1,L*K,0)
						ysub=panelsubmatrix(yhat,i,info3)
						zsub=panelsubmatrix(zhat,i,info3)
						wsub=panelsubmatrix(wvar,i,info3)
						for (j=1; j<=rows(ysub); j++) {
							yz=yz+(ysub[j,.]#zsub[j,.])*wsub[j,.]
						}
						shat3=shat3+quadcross(yz,yz)
					}
				}
			}
		}

// 1st level of clustering, no kernel-robust
// Entered unless 1-level clustering and kernel-robust
		if (!((kernel~="") & (clustvarname2==""))) {
			shat=J(L*K,L*K,0)
			for (i=1; i<=N_clust; i++) {
				yz=J(1,L*K,0)
				ysub=panelsubmatrix(yhat,i,info)
				zsub=panelsubmatrix(zhat,i,info)
				wsub=panelsubmatrix(wvar,i,info)
				for (j=1; j<=rows(ysub); j++) {
					yz=yz+(ysub[j,.]#zsub[j,.])*wsub[j,.]
				}
				shat=shat+quadcross(yz,yz)
			}
		}

// 2-level clustering, no kernel-robust
		if ((clustvarname2~="") & (kernel=="")) {
			imax=max(clustvar2)
			shat2=J(L*K,L*K,0)
			for (i=1; i<=imax; i++) {
				yz=J(1,L*K,0)
				svar=(clustvar2:==i)
				ysub=select(yhat,svar)
				zsub=select(zhat,svar)
				wsub=select(wvar,svar)
				for (j=1; j<=rows(ysub); j++) {
					yz=yz+(ysub[j,.]#zsub[j,.])*wsub[j,.]
				}
				shat2=shat2+quadcross(yz,yz)
			}
		}

// 1st level of cluster, kernel-robust OR
// 2-level clustering, kernel-robust and time is 2nd cluster variable
		if (kernel~="") {
			shat2=J(L*K,L*K,0)
// First, standard cluster-robust, i.e., no lags.
			i=min(t)
			while (i<=max(t)) {  // loop through all ts
				yz=J(1,L*K,0)
				svar=(t:==i)
				ysub=select(yhat,svar)
				zsub=select(zhat,svar)
				wsub=select(wvar,svar)
				for (j=1; j<=rows(ysub); j++) {  // sum over all obs in this t
					yz=yz+(ysub[j,.]#zsub[j,.])*wsub[j,.]
				}
				shat2=shat2+quadcross(yz,yz)
				i=i+tdelta
			} // end i loop through all ts

// Spectral windows require looping through all T-1 autocovariances
			if (window=="spectral") {
				TAU=T/tdelta-1
			}
			else {
				TAU=bw
			}
			tnow=st_data(., tindexname)
			for (tau=1; tau<=TAU; tau++) {
				ghat=J(L*K,L*K,0)
				kw = m_calckw(tau, bw, kernel)
				if (kw~=0) {						// zero weight possible with some kernels
													// save an unnecessary loop if kw=0
													// remember, kw<0 possible with some kernels!
					lstau = "L"+strofreal(tau)
					tlag=st_data(., lstau+"."+tindexname)
					tmatrix = tnow, tlag
					svar=(tnow:<.):*(tlag:<.)		// multiply column vectors of 1s and 0s
					tmatrix=select(tmatrix,svar)	// to get intersection, and replace tmatrix
													// if no lags exist, tmatrix has zero rows
					if (rows(tmatrix)>0) {
						i=min(t)
						while (i<=max(t)) {
							yz=J(1,L*K,0)
							yz1=J(1,L*K,0)
							for (j=1; j<=rows(tmatrix); j++) {
								if (t[tmatrix[j,1]]==i) {
									wvarj =wvar[tmatrix[j,1]]
									wvarj1=wvar[tmatrix[j,2]]
									yj    =yhat[tmatrix[j,1],.]
									yj1   =yhat[tmatrix[j,2],.]
									zj    =zhat[tmatrix[j,1],.]
									zj1   =zhat[tmatrix[j,2],.]
									yz=yz+(yj#zj)*wvarj
									yz1=yz1+(yj1#zj1):*wvarj1
								}
							}
							ghat=ghat+quadcross(yz,yz1)
							i=i+tdelta
						}
						shat2=shat2+kw*(ghat+ghat')
					}
				} // end i loop through all possible ts
			}  // end tau loop
// If 1-level clustering, shat2 just calculated above is shat
			if (clustvarname2=="") {
				shat=shat2
			}
		}

// 2-level clustering, completion
// Cameron-Gelbach-Miller/Thompson method:
// Add 2 cluster variance matrices and subtract 3rd
		if (clustvarname2~="") {
			shat = shat+shat2-shat3
		}		
	}  // end cluster block

	if (sw~="") {
// Stock-Watson adjustment.  Calculate Bhat in their equation (6).  Also need T=panel length.
// They define for balanced panels.  Since T is not constant for unbalanced panels, need
// to incorporate panel-varying 1/T, 1/(T-1) and 1/(T-2) as weights in summation.


		st_view(ivar, ., st_tsrevar(ivarname), touse)
		info_ivar = panelsetup(ivar, 1)

		shat=J(L*K,L*K,0)
		bhat=J(L*K,L*K,0)
		N_panels=0
		for (i=1; i<=rows(info_ivar); i++) {
			ysub=panelsubmatrix(yhat,i,info_ivar)
			zsub=panelsubmatrix(zhat,i,info_ivar)
			wsub=panelsubmatrix(wvar,i,info_ivar)
			Tsub=rows(ysub)
			if (Tsub>2) {			// SW cov estimator defined only for T>2
				N_panels=N_panels+1
				sigmahatsub=J(K,K,0)
				zzsub=J(L,L,0)
				shatsub=J(L*K,L*K,0)
				for (j=1; j<=rows(ysub); j++) {
					yzi=ysub[j,.]#zsub[j,.]
					shatsub=shatsub+quadcross(yzi,yzi)*wsub[j]
					sigmahatsub=sigmahatsub + quadcross(ysub[j,.],ysub[j,.])*wsub[j]
					zzsub=zzsub+quadcross(zsub[j,.],zsub[j,.])*wsub[j]
				} // end loop through j obs of panel i
				shat=shat + shatsub*(Tsub-1)/(Tsub-2)
				bhat=bhat + zzsub/Tsub#sigmahatsub/(Tsub-1)/(Tsub-2)
			}
		} // end loop through i panels
// Note that Stock-Watson incorporate an N-n-k degrees of freedom correction in their eqn 4
// for what we call shat.  We use only an N-n degrees of freedom correction, i.e., we ignore
// the k regressors.  This is because this is an estimate of S, the VCV of orthogonality conditions,
// independently of its use to obtain an estimate of the variance of beta.  Makes no diff aysmptotically.
// Ignore dofminus correction since this is explicitly handled here.
// Use number of valid panels in denominator (SW cov estimator defined only for panels with T>2).
		shat=shat/(N-N_panels)
		bhat=bhat/N_panels
		shat=shat-bhat
		shat=shat*N
	} // end Stock-Watson block

	_makesymmetric(shat)

// In special cases, shat may not be psd. This uses spectral decomposition to address it.
// Follows recommendation in Stock-Watson 2008 Econometrica, Remark 8.
	if (spsd~="") {
		symeigensystem(shat,Evec,Eval)
		Eval = abs(Eval)
		shat = Evec*diag(Eval)*Evec'
	}

// Finally, calcluate vhat	
	if ((LMWald=="LM") & (kernel=="") & (robust=="") & (clustvarname=="")) {
// Homoskedastic, iid LM case means vcv is identity matrix
// Generates canonical correlation stats.  Default.
		vhat=I(L*K,L*K)/N
	}
	else {
		vhat=(iryhat'#irzhat')*shat*(iryhat'#irzhat')'
		vhat=(vhat+vhat')/2
	}

// ready to start collecting test stats
	if (allrank~="") {
		firstrank=1
		lastrank=min((K,L))
	}
	else if (nullrank~="") {
		firstrank=1
		lastrank=1
	}
	else if (fullrank~="") {
		firstrank=min((K,L))
		lastrank=min((K,L))
	}
	else {
// should never reach this point
printf("ranktest error\n")
		exit
	}

	rkmatrix=J(lastrank-firstrank+1,6,.)
	for (i=firstrank; i<=lastrank; i++) {

		if (i>1) {
			u12=ut[(1::i-1),(i..L)]
			v12=vt[(1::i-1),(i..K)]
		}
		u22=ut[(i::L),(i..L)]
		v22=vt[(i::K),(i..K)]
		
		symeigensystem(u22*u22', evec, eval)
		u22v=evec
		u22d=diag(eval)
		u22h=u22v*(u22d:^0.5)*u22v'

		symeigensystem(v22*v22', evec, eval)
		v22v=evec
		v22d=diag(eval)
		v22h=v22v*(v22d:^0.5)*v22v'

		if (i>1) {
			aq=(u12 \ u22)*luinv(u22)*u22h
			bq=v22h*luinv(v22')*(v12 \ v22)'
		}
		else {
			aq=u22*luinv(u22)*u22h
			bq=v22h*luinv(v22')*v22'
		}

// lab is lambda_q in paper (eqn below equation 21, p. 104)
// vlab is omega_q in paper (eqn 19 in paper, p. 104)
		lab=(bq#aq')*vecth
		vlab=(bq#aq')*vhat*(bq#aq')'

// Symmetrize if numerical inaccuracy means it isn't
		vlab=(vlab+vlab')/2
		vlabinv=invsym(vlab)
// rk stat Assumption 2: vlab (omega_q in paper) is nonsingular.  Detected by a zero on the diagonal,
// since when returning a generalized inverse, Stata/Mata choose the generalized inverse that
// sets entire column(s)/row(s) to zeros.
// Save df and rank even if test stat not available.
		df=(L-i+1)*(K-i+1)
		rkmatrix[i-firstrank+1,2]=df
		rkmatrix[i-firstrank+1,4]=i-1
		if (diag0cnt(vlabinv)>0) {
printf("\n{text:Warning: covariance matrix omega_%f}", i-1)
printf("{text: not full rank; test of rank %f}", i-1)
printf("{text: unavailable}\n")
		}
// Note not multiplying by N - already incorporated in vhat.
		else {
			rk=lab'*vlabinv*lab
			pvalue=chi2tail(df, rk)
			rkmatrix[i-firstrank+1,1]=rk
			rkmatrix[i-firstrank+1,3]=pvalue
		}
// end of test loop
	}

// insert squared (=eigenvalues if canon corr) and unsquared canon correlations
	for (i=firstrank; i<=lastrank; i++) {
		rkmatrix[i-firstrank+1,6]=cc[i-firstrank+1,1]
		rkmatrix[i-firstrank+1,5]=ev[i-firstrank+1,1]
	}
	st_matrix("r(rkmatrix)", rkmatrix)
	st_matrix("r(ccorr)", cc')
	st_matrix("r(eval)",ev')
// Save V matrix as in paper, without factor of 1/N
	vhat=N*vhat
	st_matrix("r(V)", vhat)
// Save S matrix as in ivreg2, with factor of 1/N
	shat=shat/N
	st_matrix("r(S)", shat)
	st_numscalar("r(N)", N)
	if (clustvarname~="") {
		st_numscalar("r(N_clust)", N_clust)
	}
	if (clustvarname2~="") {
		st_numscalar("r(N_clust2)", N_clust2)
	}
// end of program
}

real scalar m_calckw(	real scalar tau,
						real scalar bw,
						string scalar kernel) 
{
	karg = tau / bw
	if (kernel=="Truncated") {
		kw=1
	}
	if (kernel=="Bartlett") {
		kw=(1-karg)
	}
	if (kernel=="Parzen") {
		if (karg <= 0.5) {
			kw = 1-6*karg^2+6*karg^3
		}
		else {
			kw = 2*(1-karg)^3
		}
	}
	if (kernel=="Tukey-Hanning") {
		kw=0.5+0.5*cos(pi()*karg)
	}
	if (kernel=="Tukey-Hamming") {
		kw=0.54+0.46*cos(pi()*karg)
	}
	if (kernel=="Tent") {
		kw=2*(1-cos(tau*karg)) / (karg^2)
	}
	if (kernel=="Daniell") {
		kw=sin(pi()*karg) / (pi()*karg)
	}
	if (kernel=="Quadratic spectral") {
		kw=25/(12*pi()^2*karg^2) /*
			*/ * ( sin(6*pi()*karg/5)/(6*pi()*karg/5) /*
			*/     - cos(6*pi()*karg/5) )
	}
	return(kw)
}  // end kw


end

* Version notes
* 1.0.00  First distributed version
* 1.0.01  With iweights, rkstat truncates N to mimic official Stata treatment of noninteger iweights
*         Added warning if shat/vhat/vlab not of full rank.
* 1.0.02  Added NULLrank option
*         Added eq names to saved V and S matrices
* 1.0.03  Added error catching for collinearities between varlists
*         Not saving S matrix; V matrix now as in paper (without 1/N factor)
*         Statistic, p-value etc set to missing if vcv not of full rank (Assumpt 2 in paper fails)
* 1.0.04  Fixed touse bug - was treating missings as touse-able
*         Change some cross-products in robust loops to quadcross
* 1.0.05  Fixed bug with col/row names and ts operators.  Added eval to saved matrices.
* 1.1.00  First ssc-ideas version.  Added version 9.2 prior to Mata compiled section.
* 1.1.01  Allow non-integer bandwidth
* 1.1.02  Changed calc of yhat, zhat and pihat to avoid needlessly large intermediate matrices
*         and to use more accurate qrsolve instead of inverted X'X.
* 1.1.03  Fixed touse bug that didn't catch missing cluster variable
*         Fixed cluster bug - data needed to be sorted by cluster for Mata panel functions to work properly
* 1.2.00  Changed reporting so that gaps between panels are not reported as such.
*         Added support for tdelta in tsset data.
*         Changed tvar and ivar setup so that data must be tsset or xtset.
*         Removed unnecessary loops through panel data with spectral kernels
*         shat vcv now also saved.
*         Added support for Thompson/Cameron-Gelbach-Miller 2-level cluster-robust vcvv
*         Added support for Stock-Watson vcv - but requires data to have FEs partialled out, & doesn't support fweights
*         Removed mimicking of Stata mistake of truncated N with iweights to nearest integer
*         Fixed small bug with quadratic kernel (wasn't using negative weights)
*         Optimised code dealing with time-series data
* 1.2.01  Fixed bug that always used Stock-Watson spectral decomp to create invertible shat
*         instead of only when (undocumented) spsd option is called.
* 1.2.02  Fixed bug that did not allow string cluster variables
* 1.2.03  Fixed bug in code for cluster+kernel robust (typo in imported code from ivreg2=>crash)
